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Abstract. One of the key aspects governing the mechanical performance of composite materi- 
als is debonding: the local separation of reinforcing constituents from matrix when the interfacial 
strength is exceeded. In this contribution, two strategies to estimate the overall response of par- 
ticulate composites with rigid-brittle interfaces are investigated. The first approach is based on 
a detailed numerical representation of a composite microstructure. The resulting problem is dis- 
cretized using the Finite Element Tearing and Interconnecting method, which, apart from compu- 
tational efficiency, allows for an accurate representation of interfacial tractions as well as mutual 
inter-phase contact conditions. The candidate solver employs the assumption of uniform fields 
within the composite estimated using the Mori-Tanaka method. A set of representative numerical 
examples is presented to assess the added value of the detailed numerical model over the simplified 
micromechanics approach. 
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1. Introduction 

Particle- reinforced composites, and the fibrous composites in particular, present a progressive 
class of materials with a steadily increasing importance in virtually all areas of structural engi- 
neering, e.g., [7, 30]. It is now being generally accepted that one of the key factors governing the 
mechanical performance of composite materials is the interfacial debonding: the partial separa- 
tion of reinforcements from matrix phase when the surface tractions locally exceed the interfacial 
strength. Therefore, a considerable amount of research effort has been invested into the develop- 
ment of a realistic model for debonding-induced damage processes in heterogeneous media in the 
field of materials science and engineering. Two major classes of models are currently available: 
(i) computational homogenization methods and (ii) micromechanics-based approaches. In spite of 
significant advances in understanding the debonding phenomenon accomplished in recent decades, 
both modelling approaches still suffer from certain limitations. 

The essential goal of the computational homogenization is to determine the detailed distribution 
of local fields within a characteristic heterogeneity pattern of the analyzed composites, usually 
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represented by a Periodic Unit Cell (PUC), employing a suitable discretization technique. In 
the context of the debonding behavior, the most frequent and versatile approach is based on 
the primal variant of the Finite Element Method (FEM), e.g., [8, 25, 45, and references therein] 
or a mixed discretization approach such as the Voronoi Cell FEM [24]. Alternative numerical 
schemes include the Boundary Element Method (BEM)-based homogenization, e.g., [6] or coupled 
BEM/FEM simulations [15]. The interfacial behavior is, in all the previously mentioned cases, 
modeled in the form of a traction-separation constitutive law, with the interfacial stiffness as the 
basic material parameter. Such concept, however, inevitably leads to problems for the perfect 
particle/matrix bonding, which formally corresponds to the infinite value of interfacial stiffness. 
In the actual implementation, of course, the perfect bonding case is approximated using a large 
penalty-like term, deteriorating the conditioning of the numerical problem manifested in spurious 
traction oscillations and hence inaccurate prediction of damage initiation, see e.g. [2, 37] and 
references therein for additional discussion. The finite interfacial stiffness, on the other hand, 
allows for mutual interpenetration of individual constituents, resulting in a non-physical distribution 
of some of the mechanical fields within composite. To circumvent these problems, an Uzawa- 
type algorithm combined with the FEM or BEM discretization were employed by Prochazka and 
Sejnoha in [33, 34] to solve a single fiber pull-out problem including interfacial Coloumb friction. 
Moreover, the primal approaches of computational contact mechanics were used with some success 
in [40, 46] to homogenize debonding composites. All these studies, however, suffer from an increased 
computational cost due to non-optimal convergence rates and as such are not well-suited for the 
fully coupled FE^ computational homogenization [14, 21], where the detailed PUC response serves 
as an "effective" non-linear constitutive relation seen at the coarser scale of resolution. 

The micromechanics-based solvers are, on the other hand, typically based on a specific solution of 
field equations and limited information of a microstructure, such as volume fractions of individual 
phases. This makes the analytical approaches extremely efficient from the computational point of 
view; their accuracy is, however, to a major extent influenced by validity of the adopted assump- 
tions. Early contributions to the micromechanical modelling of debonding composites include the 
contribution of Benveniste [3] and Pagano and Tandon [29] , where the original uniform field theo- 
ries were extended to treat possible field discontinuities at the fiber-matrix interphase. Additional 
improvements include introduction of compliant interfaces with linear behavior [18] , later extended 
for non-linear constitutive relations [9, 38]. In addition, layered inclusion models combined with 
simple numerical procedures were introduced in [39] to simulate specific particle-reinforced com- 
posites. Several attempts have also been made to increase the predictive capabilities of effective 
media models by calibrating their response against the results of detailed numerical simulations. 
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e.g. [19, 35, 42]. However, lacking reliable reference numerical solution, due to the reasons summa- 
rized above, such a strategy can hardly reproduce the dominant features of the complex interaction 
mechanisms in the material systems under investigation. 

In this contribution, we present an alternative numerical approach to the numerical homogeniza- 
tion of debonding composites, which circumvents the previously mentioned obstacles. The method 
itself builds on the Finite Element Tearing and Interconnecting (FETI) solver due to Farhat and 
Roux [13] and adopts recent ideas in duality-based solvers [11, 23, 41] to the computational homog- 
enization setting. The added value of the duality- and FETI-based framework can be summarized 
as follows: 

• For the perfect complete bonding case, the method reduces to the original FETI algorithm; 
no artificial penalty-like stiffness is therefore needed to enforce the displacement continuity. 

• Since interfacial tractions are introduced independently from the displacement fields in the 
form of Lagrange multipliers, they are captured rather accurately in the whole loading 
regime. Therefore, the interfacial strength criteria can be directly adopted, cf. [16]. 

• In addition, the Lagrange multipliers allow us to efficiently treat frictionless contact condi- 
tions [10], making the numerical algorithm well-suited for providing reference solutions to 
micromechanics schemes. 

• As typical of FETI-based methods, the dual problem involves the interface- related quantities 
only. This not only reduces the problem complexity by an order of magnitude, but also 
introduces a natural coarse level to the numerical algorithm [10, 12], which opens the way 
to efficient and scalable iterative solvers. 

In the rest of this work, these claims are clarified in more details by addressing the fundamen- 
tal case of composites with rigid-perfectly brittle interfaces in the two-dimensional setting. The 
paper starts with a brief overview of the first-order periodic homogenization of composites with 
imperfect bonding of phases in Section 2. The numerical resolution of the unit cell problem using 
the FETI method is covered in Section 3, while a simple micromechanical model based on the 
Mori-Tanaka method is introduced in Section 4. Performance of both approaches is compared 
in Section 5 for both regular as well as disordered real- world microstructures. Finally, the most 
important results are summarized in Section 6. 

In the whole text, representation of the symmetric second and fourth-order tensors of the Mandel 
type is systematically employed. In particular, a, a and a denote a scalar value, a vector or a matrix. 
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respectively and the matrix representations are defined, in the two-dimensional setting, via 
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cf. [27, Section 2.3]. Other symbols and abbreviations are introduced in the text as needed. 



2. Overview of periodic homogenization 

In the current Section, the basic notation and essentials of the first-order periodic homogenization 
theory are briefly summarized, with the aim to provide a continuous setting compatible with the 
FETI-based discretization discussed next. For more detailed treatment of individual steps, an 
interested reader is referred to [17, 26]. 




Figure 1. Decomposition of a PUC; (a) Unit cell, (b) matrix phase, (c) inclusion (fiber). 



2.1. Geometry of PUC and averaging. Consider a Periodic Unit Cell (PUC) representing a 
cross-section of composite material with long fibers, see Figure 1(a). Formally, the PUC is un- 
derstood as an open set n C with boundary F. Within the PUC, we distinguish n disjoint 
sub-domains il^*). In the following text, the index i = 1 will be reserved for the multiply-connected 
matrix phase, whereas simply-connected non-overlapping heterogeneities, i.e. fibers, are enumer- 
ated by i = 2, 3, . . . , n. It will also prove useful to introduce the internal interfaces (denoted by 
square brackets instead of round ones): 

n 

pb'] =W)n Fill = y F^] , (2) 

i.e. the boundary P'-'] corresponds to the interface between j-th fiber and the matrix phase, while 
pt^l is reserved for the whole internal matrix interface, cf. Figures 1(b) and l(c).^ The normal to 

"'^Note that in order to unify the notation, index i consistently ranges from 1 to n, while j € {2, 3, . . . , n}. 
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the set J7, defined on the whole boundary T, is referred to as n with the appropriate adjustment 
for internal boundaries. Note that for any y G the associated normals for the matrix phase 
and the j-fiber verify 

n[il(y) = -nb1(y), (3) 

cf. Figures 1(b) and 1(c). 

Finally, we fix a nomenclature related to the determination of averages for possibly discontinuous 
fields defined on PUC. To that end, consider a collection of functions /^*^ defined independently on 
sub-domains jointly denoted as 

' (y) , ye f^(') 



f{y) 



(4) 



^ /(") (y), ye 

Then, given two functions / and g, we introduce the averaging relation in the form, cf. [26, Sec- 
tion 2] : 

/df{y) , A I f 5/W(y) 



\ dy, 



•k 



where stands for the area of PUC and, in analogy with Equation (2), and denote the 
traces of functions /(*^ and g^'^^ on F^, e.g. [36]. Note that for rj-continuous functions / and g, the 
second term of the right hand side of (5) vanishes due to identity (3), thus recovering the standard 
format of the average on the unit cell. 

2.2. First-order homogenization framework. In the current work, the strain controlled ap- 
proach to the first-order computational homogenization (in the terminology introduced in [26]) is 
adopted. Within this setting, the analysis is split into two independent levels: the microscale cor- 
responding to the behavior on the level of individual constituents and the macroscale providing the 
overall response of the heterogeneous material under investigation. In the sequel, we restrict our 
attention to the microscale problem and consider a PUC subject to a macroscopic strain which, 
when combined with the governing equations of continuum mechanics and the constitutive descrip- 
tion of individual phases, determines the distribution of microscopic fields within the PUC. The 
average value of the microscopic stress cr then yields the value of macroscopic stress XI, implicitly 
defining the homogenized constitutive law. 

Following this conceptual lead, we introduce an additive split of displacement and strain fields 
u and e in the form: 

u^'\y) = MvVe + u*^'\y), £«(?/) = 5««(y) =E + e*^'\y), (6) 
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where the first part of the decompositions (6) corresponds to an affine displacement field due to 
E, whereas the second part is the correction due to the heterogeneity of the PUC. The matrices X 
and d are defined in slightly different form than in the standard approach [5] 
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due to the adopted Mandel notation. 

In general, two conditions need to be satisfied to arrive at a thermodynamically consistent macro- 
micro coupling. We start from the macro-micro strain compatibility condition 

{<y)) = E, (8) 

transformed using Equation (6) and the integration formula (5) into a generalized boundary con- 
ditions 

{e*{y))=Q^ j^u{y)u*{y)dT = Q, (9) 

imposed on the fluctuating displacement field u* . The v matrix appearing in (9) stores the com- 
ponents of the normal vector and is defined analogously to (7)i: 

ni{y) 

(10) 



V2 



n2[y) 



n2{y) 
j-,n,{y) 



The second ingredient of the macro-micro scale transition is provided by the Hill lemma 

i;Ti] = (£(y)Mi/)), (11) 

imposing the equality of the work of macroscopic and microscopic stresses and strains. Assuming 
a self-equilibrated microscopic stress field (i.e. d^a = 0) and employing the identity (5), the Hill 
lemma can be transformed into 

(s*{y)'^a{y)) = ^ jf u*{y)''X{y) dT = 0, (12) 

where A denotes the traction vector defined via 

A(y) = ^(y)M?/)- (13) 

Note that the two conditions (9) and (12) represent the necessary constraints for a consistent 
macro-micro scale tying. A particularly convenient choice, especially for the adopted periodic 
setting, is provided by restricting the heterogeneous displacements u* to Jl-periodic fields. Then, 
for any homologous vectors y^ and y~ located on the boundary T, cf. Figure 1(a), the periodicity 
and anti-periodicity of displacement and tractions holds: 

u*{y+)=u*iy-), Xiy+) = -Xiy-), (14) 
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which, together with the identity n{y^) = —n{y^), ensures the satisfaction of both kinematic and 
energetic consistency conditions. 



2.3. Unit cell problem. The distribution of the displacements and interfacial tractions fohows 
from the total energy functional in the form: 



1=1 



(15) 



where the symbol L*-*^ stands for the domain-wise constant material stiffness matrix, a is used to 
denote the trial value of a quantity a; u^^'^ and \^^\y) correspond to displacement and traction 
fields expressed in the local coordinate system as 

uf{y) 



u^^iy) 






4\y) 




_ ufiiy) 


al^iy) 
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TH(y) 
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(16) 
(17) 



XH( 



y) 



(18) 



with T[*1 denoting the transformation matrix related to the i-th interface: 

nW(y)"^ 
tW(y)T 

cf. Figures 1 and 2(b) for an illustration. Employing the displacement and strain decomposition (6) 
together with the interfacial equilibrium conditions 

AW(y) = -Ayi(?/) VyerlJl (19) 

leads to a modified energy functional 

e (u*{y),X^'\y)) = jz{\l e*«(y)TL«r«(y)dl7 + f F «(y)TLW£; dl^ 



1=1 

n 



i=2 



riji 



X^\v)'^ (u*^^Hy)-uf\y)) dr. 



(20) 



defined on the set of kinematically admissible displacements 

K = {u*{y) is 0-periodic, u*{y) = G To} , (21) 

with Fd denoting a part of the boundary, typically the corner nodes of a PUC, where the Dirichlet 
boundary conditions are imposed to prevent the rigid body modes. The admissible interfacial 
tractions are constrained to the statically admissible set A discussed in more detail in Section 3.2. 
The "true" displacement and traction fields then coincide with the saddle point of the functional: 

thereby defining the unit cell problem for composites with debonding phases. 



( u* {y) , X^^\y)) = arg max min Q (u*{y),X^\y) 



(22) 
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3. FETI-BASED SOLUTION TO UNIT CELL PROBLEM 



This Section presents an engineering approach to the solution of non-hnear debonding prob- 
lem (22). In particular, a two-step iterative procedure is employed, where in the outer loop, the 
static admissibility constraint a' ^ € A is dropped to arrive at the prediction of traction distribution 
a'"*^^. Next, in the correction step, the trial tractions are transformed into a physically admissible 
value using an interfacial constitutive law and the loop is repeated until convergence is achieved. It 
is fair to mention that, when compared to more sophisticated techniques based on quadratic pro- 
gramming [10, 11, 41], the current technique is definitely less robust due to the missing convergence 
theory. The definite advantage of the algorithm, on the other hand, is its simplicity; as illustrated in 
the sequel, the implementation systematically employs a standard FETI solver for linear elasticity, 
see e.g. [22, Chapter 4]. In the sequel, individual steps of the algorithm are summarized in some 
detail. 



3.1. Prediction step. Following the basic idea of the FETI-based discretization [12, 22], the true 
fields of fluctuating displacements and the corresponding strains are approximated independently 
on individual domains rj'-*-' 



where N„ denotes the matrix of basis functions, B„ = SN^ is the displacement-to-strain matrix 
and du stores the nodal displacements of individual components, e.g. [5]. The discretization of the 
interfacial tractions in the predictor step is performed analogously: 



Following the standard Ritz-Galerkin procedure, identical basis functions are used for the dis- 
cretization of the trial values. 

After introducing the approximations (23) and (24) into the energy functional (20), the station- 
arity conditions for du'^ and dx'^^ reduce to a system of linear equations 




Vy G 



(23) 




(24) 




fii) _ £(^^dx 



(25) 
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(26) 



1=1 
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with the individual matrices expressed in the form: 

= / B«(y)TL«B»(y)dO, (27) 

= - f £;'^L«B«(y)dl^, (28) 

n „ 

=-Y. TW(y)TNW(y)TTW(y)N«(y)dr, (29) 

gU) = I TW(y)TNW(y)TTW(2/)N(f'Hz/)dr. (30) 
The solution of Equation (25) exists only if and only if 

RWT(^y.»_^»TJ^W) =0, (31) 

where R^*^ stores the rigid body modes of the i-th domain. 

Note that in the actual implementation, the consistent compatibility matrices £^^'^ appearing 
in Equations (29) and (30) are replaced with the lumped Boolean sparse matrices, which enforce 
the displacement continuity discretely at the corresponding nodes of the finite element mesh and 
furnish dx ^ with the physical meaning of the concentrated nodal forces, see [22, Section 4.2] for a 
more detailed discussion. 

Now we proceed with expressing the coefficients of fluctuating displacements du from the systems 
of equations (25) in the form 

d« = K«t (/» - + R«d» (32) 

The first term in relation (32) corresponds to the particular solution of the i-th component of the 
system (25), which is expressed using the generalized inverse matrix K^*)''^ replacing the inverse 
matrix for singular K^*). The second term then corresponds to a homogeneous solution, expressed 
as the linear combination of rigid body modes R*^*) with coefficients d^^ to be determined. Next, 
we substitute the coefficients of fluctuating displacements di*^ from relations (32) to system (26) 
and the solvability conditions (31) to account again for a potential singularity of matrices K^*): 
^^WK»t£:WTrf~^W _ Y^£{i)^{i)df = ^£:WKWty.W, (33) 

i=l 1=1 i=l 



Observe that the elimination of the primary unknowns du leads to a substantially smaller dual 
problem, formulated in terms of variables d\^'^ and d^^ . Such system can be efficiently solved using 
the Modified Conjugate Gradient (MCG) method, augmented by the projection step to enforce the 
solvability condition (31). As will be reported separately, a special structure of the unit cell problem 
can be efficiently exploited to construct the pseudoinverses K*^*^^ efficiently. The robustness and 
scalability of the numerical algorithm can further be extended by applying orthonormalization of 
matrices R^*^ and appropriate preconditioning procedure. 
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3.2. Interfacial constitutive law. To close the problem statement, an appropriate interfacial 
constitutive model needs to be specified. The particular choice adopted in the current work ne- 
glects the interaction of the normal and tangential strengths Cmax and ''"max and accommodates the 
tangential slip in presence of compressive normal stresses. Such assumptions lead to a three-state 
description appearing in Figure 2(a). In particular, the following modes are distinguished: (i) per- 
fect bond corresponding to an intact interface, (ii) frictionless contact state, where the local shear 
strength is exceeded and sliding between fiber and matrix is activated and (iii) complete separation 
of phases. 







F contact 
I problem 


complete 
separation 


f perfect 




bond 






~ '^max 




(a) (b) 
Figure 2. Interfacial constitutive law; (a) three-state model, (b) local coordinate system. 



The set of admissible tractions then specializes to 

A = < ^max n |fW(y)| < w) U < o) Vy G fW} . (35) 

Correspondingly, the constraint A E A is simply enforced by setting the traction value to zero 
whenever the corresponding strength is exceeded by the predicted value A. Note that it follows 
from duality theory that the condition a'^^ < is equivalent to the non-interpenetration for the 
matrix and fibers (un' > tin'), cf. [11]. 

3.3. Implementation strategy. The resulting incremental algorithm compatible with the 
adopted interfacial law (35), inspired by a related study [46], is briefly summarized by means 
of a pseudo-code shown in Table 1. Note that in order to simplify the exposition, a proportional 
loading path in the form 

E(t) = tE^^^ with 0<t<l (36) 

is considered. Moreover, the irreversibility of debonding is enforced using the static condensation [5, 
Section 2.5] of the corresponding entries of d}-^\ cf. lines 8 and 12 of Table 1. 
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i 


initialize material parameters (phase elastic properties) and geometry (mesn) 








assemble stiffness matrices K*^'\ 


Eqs. 


(27) 




transformation matrices T'' and compatibility matrices ' 


Eqs. 


(29,30) 


li 


Initialize loading data: macroscopic strain .c/max 








assemble load vectors j„/ax 


Eqs. 


(28) 


iii 


Computation 1:* _h 1 1- based quantities: 
pseudomverse and rigid body R'- ^ matrices 






iv 


betup mterlacial strengths and pseudo-times ti < 12 < ■ ■ ■ < tj^ 






1 


tor fc = 1 : iV (cycle through time steps) 






2 


set f-'' = ifc/,V/ax and £; = tfeiJ^ax 






3 


do (outer loop) 






4 


predict distribution of mterfacial tractions a\ using MCG 


Eqs. 


(33,31) 


5 


tor n — 1 . M (internal cycle though mterlacial nodes) 






6 


extract a and f for n-th node 






7 


it (T < U and \t\ > Tmax (contact problem) 






8 


set a — a and r = by static condensation of ' (without possibility of healing ) 






9 


elseit cr < fJmax and \t\ < Tmax (pcrlcct bond) 






10 


set a — a and r = r 






11 


else (complete separation) 






12 


set CT = and r = by static condensation of (without possibility of "healing") 






13 


endit 






14 


endfor 






15 


until interfacial conditions remain unchanged 






16 


endfoi 







Table 1 . Conceptual implementation of FETI-based solver for debonding problem 



4. MiCROMECHANICAL MODEL 

In the current Section, a simplified micromechanics-based solution to the debonding problem 
is constructed, in order to, at the same time, verify the developed FETI-based algorithm as well 
as to assess the added value of the the detailed numerical modelling. To that end, a composite 
is treated as a two-phase material, consisting of a matrix phase (domain ^2'-"'^^), occupied by (n — 
1) indistinguishable particles (related to domains O^^^, . . . , il^") in the FETI-based method). 
Moreover, the available geometrical information is reduced to the phase volume fractions c^^^ and 
c^'^\ defined as 
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where, in accord with the rest of the current section, a symbol "^^^ is reserved for a quantity related 
to the matrix phase while •'^^^ collectively describes particle-related unknowns. 

Employing the simplified data, the overall stiffness of a two-phase composite reads 

L^ff = L(i) + c(2) - L(i)) A(2), (38) 

where L^^^ and L^^-* are the matrix and particle stiffnesses, respectively, and A^^^ stands for a strain 
concentration factor, relating a macroscopic strain E to the average strain in a particle. It follows 
from the Benveniste's reformulation [4] of the original Mori-Tanaka method [28] that the strain 
concentration factor can be accurately approximated as 

A(2) « Ai^^ = T(2) (c(i)l + c(2)t(2)) , (39) 

with T^^) denoting the partial strain concentration factor (with a physical meaning of the average 
strain in a particle due to prescribed average strain in the matrix), expressed in the form 

t(2) = (i + P(2)(^L(2)-lW))"\ (40) 

The term P^^^ appearing in Equation (40) is the Walpole matrix [44], related to a response of an 
isolated inclusion embedded in an infinite matrix phase, which, under the assumption of circular 
fibers, isotropic matrix and plane strain state, can be explicitly expressed as, cf. [44]: 

8z^(i) - 5 1 

1 8i/(i)-5 , (41) 

32i/(i) - 24 

where E^^'^ a u^^^ are the Young modulus and the Poisson ratio of the matrix phase. 

When assuming a phase- wise uniform distribution of mechanical fields within heterogeneities, 
the stress value in a particle can be estimated as 

a('\y)^a^^^ = L(')A^^^E, (42) 

and, employing the equilibrium conditions (19), converted into the approximate value of boundary 
tractions at the particle/matrix interface 

aW(0) « aWmt(0) = ivW(0)Tcr(2), (43) 

where 9 is an angle parameterizing a position on the interface introduced in Figure 2(b) and 
the stress-traction transformation matrix u^^'^ is expressed using Equation (10) with n = n^^'^ = 
— QQgQ —sin 6 The algorithmic implementation of the micromechanics-based solver closely 
follows the procedure introduced by Table 1. The debonding mechanism is reduced to a two-state 
description, with one state describing an undamaged composite, whereas the second state corre- 
sponds to completely debonded particles. In particular, the perfect bonding at the matrix/particle 
interface is maintained until the normal or tangential tractions do not exceed the local strengths, 

(T^Ti^) < '^max and |rj5!p(6')| < w for ah < < 27r, (44) 



p(2) = 1 + -^'^ 

8ew [i^m - 1) 
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where cr^j, and Tj^i^ denote the normal and tangential tractions determined from (43) via rela- 
tion (17) with t[^l = sin 9 — cos 6 """.If the previous condition is violated for at least one value 
of 9 for a given level of macroscopic strain, we irreversibly set L^^^ = 0. 

5. Numerical examples 

Performance assessment of the numerical homogenization and the micromechanics-based solvers 
will be performed on the basis of two PUCs representing cross-sections typical of unidirectional 
composites: a parametric hexagonal unit cell model analyzed in Section 5.1 and a twenty-particle 
unit cell treated in Section 5.2, which was obtained in [47] by a careful analysis of a graphite/polymer 
composite sample. All FETI-based examples presented bellow are determined using an in-house 
code derived from a Matlab implementation of finite element solver introduced in [1] with the 
DistMesh code, developed by Persson and Strang [31], used to generate finite element meshes. The 
material properties of individual phases and interface are taken from a related study [19] and appear 
in Table 2. The plane strain assumption is adopted to reduce the analysis to a two-dimensional 
cross-section and the displacement fields are discretized using constant strain triangular elements 
with the characteristic element /fiber diameter ratio approximately equal to 0.05. In all reported 
numerical simulations, the following discretization of the time interval [0, 1] was considered. First, 
assuming perfectly bonded phases, the value t* related to the debonding onset was determined. 
Next, the post-peak branch (t*, 1] was uniformly sampled using N equisized intervals. An interested 
reader is referred to [17] for additional implementation-related details. 



Matrix 


Young's modulus, E^^^ 


1 GPa 


Poisson's ratio, z^^^^ 


0.4 


Inclusion 


Young's modulus, E^"^^ 


150 GPa 


Poisson's ratio, v^"^^ 


0.3 


Interface 


Normal strength, cJmax 


0.02 GPa 


Tangential strength, Tmax 


0.02 GPa 



Table 2. Material data of composite system 



5.1. Hexagonal packing of particles. As the first example, we consider a hexagonal PUG sub- 



ject to a bi-axial tensile/compressive strain load characterized by -En 



0.02 



-0.02 
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(a) (b) 
Figure 3. Scheme of debonding initiation for a hexagonal PUC; (a) N = 100, 
(b) N = 10. Displacements are scaled by factor of 20. 



recall Equation (36). Stress paths for the MT scheme correspond to a fine discretization of the 
overall time {N = 100), whereas the value = 10 was adopted for the FETI-based algorithm. 
Note that the coarse time stepping of the latter approach is mainly related to a high sensitivity 
of the post peak branch to inevitable geometrical imperfections of the finite element mesh. As 
illustrated in Figure 3, taking N = 100 leads to a non-symmetric debonding mode, resulting in a 
substantially different distribution of local fields and overall response than for the symmetric defor- 
mation mechanism. When the time discretization is set to = 10, on the other hand, the larger 
value of the overall strain increment suppresses the mesh-induced imperfections and the solution 
symmetry is maintained. 

The resulting macroscopic stress-strain relations obtained by the MT and FETI-based approaches 
for the fiber volume fraction c^^^ = 50% appear in Figure 4. For comparison, the ideal debonding 
mode, corresponding to sufficiently small time increment and perfectly symmetric mesh, is indicated 
by gray color. 




Figure 4. Macroscopic stress-strain diagram for hexagonal PUC; c^^) = 50%. 
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In the elastic regime, both methods predict almost identical values of macroscopic stresses with 
the MT solver providing a slightly more compliant response than the FEM-based approach. Such 
behavior is consistent with the fact that the Mori-Tanaka estimates, for material systems with 
stiffer particles embedded in a weaker matrix phase, coincide with the lower Hashin-Shtrikman 
variational bounds, cf. [32]. A somewhat higher discrepancy is observed for the prediction of a 
debonding onset t* , where the "unsafe" MT result yields ij^rp = 0.41, which by about 25% exceeds 
the reference value tpETi ^ 0.33. The difference can be attributed to the value of constant fiber 
stress field adopted in the MT estimate. Such an assumption is well-justified when determining 
the mean response of a material with coherent interfaces, as demonstrated by excellent match of 
the elastic branches in Figure 4, but becomes less accurate when estimating the extreme stresses 
governing the failure process. This claim is further supported by plots of interfacial tractions at 
the debonding onset shown in Figure S.'^ We observe that, for the current case of 50% particle 
volume fraction, the deviation of the fiber stress from the assumed uniformity in the MT method 
is further propagated to the values of local tractions. It is also worth noting that, up to some 
discretization-based effects, the FETI solution is free of spurious traction oscillations appearing 
when introducing interface elements between fibers and matrix as reported, e.g., in [37]. 




Angle d [deg] Angle 9 [deg] 

(a) (b) 
Figure 5. Distribution of interfactial tranctions at the onset of debonding; (a) nor- 
mal and (b) tangential components; discrete values correspond to the FETI solution, 
the MT results are independent of c*-^^ . 

After the debonding initiation, the discrepancies between both approaches become even more 

pronounced. The MT solver predicts, due to the assumption of complete interfacial failure, a 

sudden drop of the stress values in both directions, reaching the values corresponding to an isotropic 

^Note that the values of tractions shown in Figure 5 correspond to diflerent debonding times t* and therefore also 
to different loading levels. 
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porous medium in the post-peak regime. The FETI-based model, being capable of capturing partial 
interfacial debonding, leads to a highly anisotropic mechanism: in the compressive direction, the 
response remains almost insensitive to the debonding phenomena, whereas the actual stress drop 
magnitude in the yi direction exceeds the value predicted by the MT scheme, leading eventually to 
the appearance of overall compressive stresses Sn due to intra-phase contact. As these conclusions 
remain valid for both "ideal" and actual discretized stress paths, we will perform the assessment 
on the basis of the latter results in the sequel. 




Next, we investigate the influence of Tmax parameter on the overall behavior of the hexagonal unit 
cell, cf. Figure 6. It appears that, for the current test, the MT approach does not capture the change 
of failure mode when increasing the interfacia-l strength, from T^ax 

from 0.02 to 0.04 GPa. Such result 
can be again explained by inspecting traction distribution for c^"^^ = 50% shown in Figure 5. Due 
to the fact that -En = I-E22I) the MT theory predicts identical extereme values of normal and shear 
tractions, which, in turn, implies that the debonding initiation is governed solely by the minimum 
value of fTmax and Tmax- The slight difference in the extreme shear and tensile tractions predicted 
by the FETI solver, on the other hand, leads to a three distict stress-strain curves, corresponding 
to a shear initiation for Tmax = 0.01 GPa and Tmax = 0.02 GPa, while for Tmax = 0.04 GPa the 
debonding starts in the normal mode. Note that similar results are found when varying the cimax 
parameter. 

The last aspect to be analyzed remains the particle volume fractions, cf. Figure 7. For the Sn 
component, the results of the FETI-based solver predict a gradual transition of the debonding 
branch from tensile to compressive stress values with increasing particle volume fraction c^'^\ A 
closely related trend, resulting from stress redistribution due to intra-phase contact, is exhibited 
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by the S22 path, where the stress drop after debonding changes its sign with increasing c^^^. The 
micromechanics-based solver, on the other hand, leads to a symmetric response in terms of normal 
macroscopic stresses with the difference in debonding time ranging from 7% to 40%. This trend is 
in a perfect agreement with traction plots in Figure 5, demonstrating that the mismatch between 
traction distribution for the MT and FETI solvers generally increases for densely packed particles 
due to more complex interaction mechanism. In particular, for low volume fractions, the results 
of the analytical approach almost exactly duplicate the numerical data; the MT method, however, 
fails to reproduce the traction redistribution at increased particle volume fractions c^'^\ 



5.2. Real-world microstructure. The universality of the conclusions reached for the hexagonal 
packing of particles is further supported by an analogous analysis of a twenty-fiber unit cell with 
particle volume fraction c^^^ = 43.9%, subject to a pure shear strain loading path with .Emax = 
0.02\/2 The mutual comparison of macroscopic stress path appears in Figure 8, storing 
the results of the MT method and, to illustrate the effect of microstructural modeling, of the 
FETI solver for the disordered microstructure and the regular hexagonal packing with identical 
fiber volume fractions. 

The stress paths predicted by FETI reflect somewhat different failure mechanisms for disordered 
and regular composite. In particular, the inelastic path for the 20-fiber unit cell exhibits a piecewise 
linear evolution of the shear stresses, accompanied by a monotonically descreasing normal stresses. 
The observed deformation mechanism can be explained by the distribution of the local stress fields 
captured in Figure 9. In particular, the debonding onset is driven by stress amplification in the 
form of diagonal shear bands visible in Figure 9(a). After exceeding the interfacial strengths, 
see Figure 9(c), the concentrations are released at newly formed discontinuity interfaces, resulting 
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in the formation of shifted shear bands connecting non-debonded particles and in local matrix- 
fiber contacts. Such mechanism is progressively repeated until the fully debonded state is reached, 
cf. Figure 9(c). For the regular microstructure, the debonding initiates simultaneously at all fibers 
and the stiffness immediately reaches the residual value (Figs. 9(b), 9(d) and 9(d)). It is worth 
noting that even the debonded composites we observed (up to minor discretization-induced effects) 
the stress values Sn = S22, which confirms the spatial isotropy of the samples in question. 

The proposed simple MT-based solution, on the other hand, shows a more compliant residual 
stiffness in the shear mode and predicts appearance zero normal stresses. Such response is again a 
direct consequence of the "porous" approximation adopted in the post-peak regime. Finally observe 
that the debonding onset pseudo-times, as predicted by different simulations, verify 

*FETI(20) = 0.19 < iFETI(hex) = 0-32 < ^MT = 0-46, 

cf. Figure 8, eventually demonstrating that the differences between the regular and irregular mi- 
crostructures due to local stress amplicifactions can be well comparable with the value found when 
comparing the analytical approach to the detailed numerical model. 

6. Conclusions 

In the current work, an overview of the FETI-based procedure for the homogenization of debond- 
ing composites was presented. In addition, the developed algorithm was used to perform a detailed 
assessment of a simplified micromechanics-based solver based on the MT method. The results of 
the performed studies have led to the following conclusions: 

• The duality-based solver offers very promising approach to the homogenization of debond- 
ing composites, both from the point of view of numerical efficiency as well as the direct 
treatment of traction-based interfacial constitutive laws. 



0.02, 




Loading parameter t 

Figure 8. Macroscopic stress-strain response for 20-fiber PUC. 



HOMOGENIZATION OF COMPOSITES WITH INTERFACIAL DEBONDING 



19 




J f 



/ 



J r 



-0.01 -0.005 0.005 0.01 0.015 0.02 0.025 



(a) t = 0.19 (displacements are scaled (b) t = 0.32 (displacements are scaled 



by factor of 20) 



by factor of 20) 





-0.025 -0.02 -0.015 -0.01 -0.005 0.005 0.01 



(c) t = 0.27 (displacements are scaled (d) t = 0.39 (displacements are scaled 
by factor of 20) by factor of 20) 




-0.06 -0.05 -0.04 -0.03 -0.02 -0.01 0.01 



(e) t = 1 (displacements are scaled by (f) t = 1 (displacements are scaled by 



factor of 5) 



factor of 5) 



Figure 9. Distribution of maximum principal stress for a 20-fiber and equivalent 
hexagonal PUCs. The values are given in GPa. 
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• The simple micromechanics-based solver is fully capable of accurate prediction of the elastic 
part of the loading branch. The onset of debonding process is captured with the difference 
from the numerical results ranging from approximately 10% to 40% for regular packing of 
particles. 

• The assumption of completely debonded interfaces adopted in the MT solver, on the other 
hand, appear to be quite unrealistic and is mainly responsible for the inconsistent predic- 
tions in the post-peak regime. 

• The results presented in Section 5.2 clearly confirm that the apparent post-critical response 
is influenced not only by details of particle arrangement, but also by the PUC size, cf. [43]. 
When dealing with inherently random real world composites, the choice of an appropriate 
statistically equivalent PUC becomes rather delicate and generally requires a proper statis- 
tical characterization in terms of geometry and mechanical quantities of interest, see [20, 48] 
for additional discussion. 

The future extension of the FETI-based solver would involve a treatment of more realistic inter- 
facial constitutive laws. For the micromechanics-based solver, the essential feature to incorporate 
seems to be the debonding-induced anisotropy in composite including the effect of contact stresses. 
Both topics enjoy our current interest and will be reported on separately. 
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